Temperature dependent surface relaxation for Al(llO) and Mg(lOlO) studied by orbital 

free ab initio molecular dynamics. 
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We have performed orbital free ab initio molecular dynamics simulations in order to study the 
thermal behaviour of two open surfaces of solid metallic systems, namely the (110) face of fee Al and 
the (lOiO) face of hep Mg. Our results reproduce qualitatively both the experimental measurements 
and previous ab initio calculations performed with the more costly Kohn-Sham approach of Density 
Functional Theory. These calculations can be viewed as a validation test of the orbital free method 
for semiinfinite surfaces, and the results underpin its reliability. 
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I. INTRODUCTION 

If a bulk metallic crystal at zero temperature is in- 
stantaneously separated into two halves, exposing two 
pristine surfaces, then the electrons and the ions redis- 
tribute themselves, responding to the new environment, 
in order to lower the total energy. In the simplest cases 
this results in surface relaxation, where the spacing be- 
tween layers near the surface varies with respect to its 
bulk values Screening of the surface and smoothing 
of the charge usually leads to an expansion of the first 
interlayer distance for close-packed faces. Inner inter- 
layer distances usually change very little for this type of 
surfaces. For open surfaces, in addition, there is a pos- 
sibility of diminishing the undercoordination felt by the 
outermost atoms by moving closer to the second layer, 
resulting usually in a contraction of the first interlayer 
distance. Inner interlayer relaxations also occur, in gen- 
eral leading to a damped oscillatory pattern of expansion, 
contraction, expansion and so on. 

When the temperature is raised additional dynamic ef- 
fects take place and it is possible to find experimentally 
a widely varying thermal behaviour. While general theo- 
retical considerations conclude that thermal expansion is 
expected 0, some surfaces exhibit an anomalously large 
effect for its first interlayer distance ( Pb(llO), Ni(lOO), 
Ag(lll), Cu(llO), Be(OOOl) ) '^l and some others show a 
negative thermal expansion coefficient for this first inter- 
layer distance, followed by positive coefficients for inner 
interlayer distances ( Al(llO) ) 3 or even by an alter- 
nating sign behavior ( Mg(lOlO)' 5], Be(lOiO) 6]). 

In this work we will focus on two of these "anoma- 
lous" systems, namely Al(llO) and Mg(lOlO). Both 
are simple sp bonded metals, amenable to study within 
a pseudopotential formalism. Moreover, both of them 
have been studied by ab initio methods, namely using 
the Kohn-Sham (KS) Q version of Density Functional 
Theory (DFT) 8]. For Al(llO) ab initio molecular dy- 
namics (AIMD) simulations were performed 9J, which 
were able to reproduce both the oscillatory relaxations 
and the thermal behaviour observed experimentally. In 
the case of Mg(lOlO) theoretical calculations based on 



the quasiharmonic approximation (QHA) and static KS- 
DFT computations |5|], reproduced also the oscillatory 
pattern both in relaxations and in thermal expansion co- 
efficients. It should be noted, however, that no AIMD 
simulations have yet been performed for this system. 

KS-AIMD simulations represent a very powerful tool 
to study metaUic surfaces, because, due to the use of 
DFT, the response of the ions to the rapid decrease of 
the electron density near the surface is calculated selfcon- 
sistently. However its application has been very scarce 
in the literature because of practical difficulties, namely 
they are extremely expensive computationally. Some of 
this cost can be alleviated if one returns to the original 
formulation of DFT which uses the electron density 
as the only variable in the theory, without any resort to 
KS orbitals. In such an orbital free (OF) theory [loj . 
the electronic kinetic energy must be computed approx- 
imately (instead of exactly in the KS version), and lo- 
cal pseudopotentials are needed (in the KS version also 
nonlocal pseudopotentials can be used). However, the 
computing time and memory saved by disposing of the 
orbitals can then be invested in studying larger systems 
for longer times. This, while important in solid metallic 
surfaces, shows its full potential in the study of liquid 
metallic surfaces where the absence of long range 
order requires the use of large samples in order to ob- 
tain reahstic results. Note that OF- AIMD still use DFT, 
and therefore the main power of the AIMD is preserved, 
i.e. the electrons and the ions near the surface respond 
selfconsistently to the rearrangement of one another. 

While OF- AIMD simulations have been successful in 
the study of static and dynamic properties of bulk liquid 
metals [iflliaUJUj, as well as in the understanding of 
the thermal properties of some metallic clusters (15| (in- 
cluding an anomalous variation of the melting tempera- 
ture in Na clusters with size 116]), its more approximate 
character (as compared to KS-AIMD ones) might induce 
someone to wonder if the theory is at all applicable to 
metallic surfaces, or at least to question how accurate it 
is. In this work we perform OF-AIMD simulations for 
the two solid systems mentioned previously, Al(llO) and 
Mg(lOlO), for which the KS calculations can be taken 
as a benchmark. The case of Al will permit a direct 
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comparison with both experiment and KS simulations, 
while for Mg the comparison with KS calculations is less 
direct, since no KS simulations were performed; never- 
theless the comparison with the results obtained within 
the QHA will be valuable anyhow. We will show that 
the OF-AIMD method is indeed able to reproduce qual- 
itatively the experimental and KS findings, butressing 
therefore its reliability. 

In section^we outline the theoretical basis behind the 
OF-AIMD simulations. Section IITTI shows our results for 
Al(llO) and Mg(lOlO), and we discuss these results and 
obtain our conclussions in section IWI 

II. METHOD 

A simple sp bonded metal is treated as a set of N 
bare ions with valence Z , enclosed in a volume V , and 
interacting with iVo = NZ valence electrons through an 
electron-ion potential v{r). The total potential energy of 
the system can be written, within the Born-Oppenheimer 
approximation, as the sum of the direct ion-ion coulom- 
bic interaction energy and the ground state energy of the 
electronic system, Eg^ under the external potential cre- 
ated by the ions, 14xt(r, {Ri}) = Yh=i v{\r - Ri\) , 

i^dfi/}) i5 I +EM^^y<^ArARi})\. (1) 

where pg (f) is the ground state electronic density and Ri 
are the ionic positions. 

According to DFT, the ground state electronic den- 
sity, Pgir), can be obtained by minimizing the energy 
functional E[p], which can be written 

E[p{r)] = Ts [p] + Eh [p] + E^, [p] + E,^t [p] (2) 

where the terms represent, respectively, the electronic ki- 
netic energy, Tg [p] , of a non- interacting system of density 
p{r), the classical electrostatic energy (Hartree term), 

^"H^i//*'^.-ff , (3) 

the exchange-correlation energy, i?xc[p], for which we 
adopt the local density approximation, and finally 
the electron-ion interaction energy, Eoxt[p\, where the 
electron-ion potential has been characterized by a local 
ionic pseudopotential which has been constructed within 
DFT. 




effort involved in this approach for large systems is allevi- 
ated in the OF-AIMD approach by use of an explicit but 
approximate functional of the density for [p] . Proposed 
functionals consist of the von Weizsacker term, 

Twipir^] = dr |Vp(r)|Vp(r), (5) 

plus further terms chosen in order to reproduce correctly 
some exactly known limits. Here, we have used an aver- 
age density model, where Tg = Tw + Tp, 

T^-yJ drp{r^'''~''~m' (6) 
k{f) = {2k%.f j dsk{s)wp{2k°p\f - s\) 

k(f) = (Stt^)^/'^ p{'f^^ , kp is the Fermi wavevector for 
mean electron density pe = Ne/V, and wp^x) is a weight 
function chosen so that both the linear response theory 
and Thomas-Fermi limits are correctly recovered. Fur- 
ther details are given in reference 

Another key ingredient of the energy functional is the 
the local ion pseudopotential, Vps{r)^ describing the ion- 
electron interaction. For each system, the Vps{r) has 
been constructed from first principles by fitting, within 
the same Ts[p] functional, the displaced valence elec- 
tronic density induced by an ion embedded in a metallic 
medium as obtained in a KS calculation. Further details 
on the construction of the pseudopotential are given in 
reference 0|. 

Given an ionic configuration, the electronic ground 
state is obtained, the potential energy for the ions is 
evaluated and the forces acting on them found using 
the Hellmann-Feynman theorem. These are then used to 
move the ions according to Newton equations of motion 
into a new configuration, after which the whole procedure 
is repeated. 

III. RESULTS 

In order to compare our OF-AIMD data with those 
obtained by KS calculations we have used in our sim- 
ulations exactly the same setup as in the previous KS 
studies. As a result, in the case of Al(llO) we have 8 
layers of 9 atoms each plus a vacuum of 8.5 A in our sim- 
ulation cell, in which the in-plane lattice spacing is taken 
as the experimental one for each temperature considered. 
In the case of Mg(lOiO) we consider 16 layers of 4 atoms 
each and an 8.5 A vacuum, with again the experimental 
in-plane lattice spacings. 

A. Al(llO) 



In the KS approach to DFT Tg [p] is calculated exactly 
by using single particle orbitals. The huge computational 



The OF-AIMD simulations of Al(llO) have been per- 
formed for two temperatures close to those where ex- 
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FIG. 1: Results for the interlayer relaxations in percent 
100(dy(T)-d^/"^(T))/d^/"^(T), for Al(llO). Full circles: OF- 
AIMD results. Line: best linear fit to them. Hatched squares: 
uncorrected KS-AIMD results. Lozenges with lines: experi- 
mental data. 



perimental data obtained through low energy electron 
diffraction (LEED) were reported, namely T = 70 and 
310 K, and also at a higher temperature of 707 K, near 
one of the KS-AIMD simulations. 

For comparison, the KS-AIMD simulation time was 
between 5 and 10 ps depending on the temperature (for 
T = 700 K it was 6 ps) with runs performed on a vector 
supercomputer (Hitachi S3600), and the OF-AIMD sim- 
ulation times were between 8 and 16 ps after an initial 
equilibration time of 2-4 ps (8 and 4 ps respectively for 
T = 707 K), with runs performed on an Intel Centrino 
laptop (with a Pentium M processor). 

The interlayer distances obtained from our simulations 
are plotted in figure 1, together with the experimental 
data and the results of the KS-AIMD simulations. Note 
that the KS-AIMD results published in reference Q were 
corrected rigidly by the difference between a KS calcu- 
lation at T = K using the same setup and another 
one using more layers, one atom per layer, and better 
Brillouin zone sampling. In figure 1 we have plotted the 
uncorrected results in order to make a fair direct compar- 
ison between the OF-AIMD and the KS-AIMD results for 
the simulation setup used. 



B. Mg(lOlO) 

The OF-AIMD simulations for this system have been 
performed near the temperatures at which the LEED 
study was performed, which also coincide with the tem- 
peratures of the KS-QHA calculations, namely T — 
106, 308 and 399 K. The equilibration and production 
runs spanned a time of 4 and 8 ps respectively. 



FIG. 2: Results for the interlayer relaxations in percent 
100(dy(T) - <i°f^{T))/d\f'{T), for Mg(lOlO). Full circles: 
OF-AIMD results. Thick continuous lines: best linear fit to 
them. Hatched squares: KS-QHA results. Dashed lines: best 
linear fit to them. Lozenges with error bars: experimental 
data. Thin continuous lines: best linear fit to them. 

Figure 2 shows the results for this system. For the 
hep structure in this orientation there are two types of 
interlayer distances, a short interlayer distance, between 
the first and second layers, between the third and fourth 
layers and so on, and a long interlayer distance (twice 
as large in the bulk solid) between the second and third 
layers, between the fourth and fifth layers and so on. We 
have separated figure 2 into two panels in order to appre- 
ciate more clearly the thermal variation of the interlayer 
distances of both types. 

IV. DISCUSSION AND CONCLUSIONS 

Starting with Al(llO) we first remark that the OF- 
AIMD simulations recover the oscillatory pattern of in- 
terlayer relaxation observed experimentally, and also the 
thermal variation of these data, showing a negative first 
interlayer thermal expansion coefficient and a positive 
one for all the inner interlayer relaxations. The agree- 
ment with the experimental data at 70 and 310 K is ex- 
cellent, but we would rather emphasize the reproduction 
of the trends rather than the numbers. One of the main 
reasons for stressing this is the lack of a detailed analy- 
sis of size effects in the simulations. According to recent 
all-electron first principles calculations of the properties 
of Al surfaces [13 , it might be neccessary to include as 
many as 23 layers in the simulation for the (110) orien- 
tation in order to obtain fully converged results. This of 
course needs to be tested for the OF-AIMD simulations 
before a comparison with experimental data can be made 
at a quantitative level. 

Comparing with KS-AIMD results, we find a reason- 
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ably good agreement, taking into account the differences 
in the simulations (kinetic energy functional and pseu- 
dopotentials). In any case, again we think that the im- 
portant point is the reproduction of the trends observed 
in the KS simulations. 

Coming to Mg(lOlO), we again remark first that the 
OF-AIMD results reproduce the trends of relaxation and 
thermal variation in this system. The short interlayer 
distances both contract and show negative thermal ex- 
pansion coefficient, while the long interlayer distances 
both expand and show positive thermal expansion coef- 
ficient. The magnitude of the contractions is reproduced 
with better accuracy than that of the expansions, but 
nevertheless the trend is the correct one. 

When comparing with the KS-QHA data we outline 
three points. First, the magnitude of the relaxations ob- 
tained from the KS-QHA is closer to the experimental 
one than that of the OF-AIMD simulations. Second, the 
thermal variation of the KS-QHA and OF-AIMD results 
is very similar, as observed from the slope of the linear 
fit to both types of data. And third, both approaches 
underestimate largely the thermal variation found in the 



experimental measurements, which show a much larger 
slope for the first three interlayer distances. 

Summarizing, the OF-AIMD results for the structure 
of the open surfaces considered and their thermal varia- 
tion reproduce qualitatively the experimental trends. In 
many cases, the results are also very similar to those 
obtained through more; expensive methods as KS-DFT, 
either used within AIMD simulations or within the QHA, 
with the only exception of the magnitude of the expan- 
sion of the long interlayer distances in the hep structure. 
A quantitative comparison with experimental data is at 
present not sensible, since a detailed study of size effects 
in the simulations is necessary, most surely for Al(llO) 
but probably also in the case of Mg(lOiO). This analysis 
is under way and will be presented elsewhere. 

In our opinion, the results shown in this work, added 
to those already published for bulk metallic liquids and 
for metallic clusters, further demonstrate that the OF- 
AIMD method is not only practical, but also reliable for 
the study of metallic systems, and in particular metallic 
surfaces. 
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